(* Mathematica Package *)

BeginPackage["GenCluster`", {"MathMR`"}]
(* Exported symbols added here with SymbolName::usage *)  

$WorkDir::usage="working directory";
$NBATH::usage="Number of bath spin";
$BathCoordinates::usage="coordinates of the bath";
$Threshold::usage="distance threshold.";
$LinkMap::usage="link map";

$BathIndexList::usage="";
$ConnectMap::usage="";

FullClusters::usage="";
SubIndexList::usage="";


ReadCoord::usage="Read Coordinates from a file.";
GenerateLinkMap::usage="Generate link map";

GenerateConnectionMap::usage="";
ClusterGrowth::usage="Grow new clusters from old one.";

SetLink::usage="set link value,true or false, to each link table";
IsConnect::usaage="return the bool value of linking";
FormLinkTable::usage="Generate all possible link table."
GenerateFullClusters::usage="";
GenSubIdx::usage="";


Begin["`Private`"] (* Begin Private Context *) 


ReadCoord[filename_]:=Module[{res},
$BathCoordinates=Import[$WorkDir<>filename];
$NBATH=$BathCoordinates//Length;
res=$NBATH;
Return[res];
]

GenerateLinkMap[]:=Module[{res,dist,halfmat},
halfmat=Table[
	If[i<j,
	dist=Norm[$BathCoordinates[[i]]-$BathCoordinates[[j]]];
	If[dist<$Threshold,1,0]
	,
	0
	]
	,{i,1,$NBATH},{j,1,$NBATH}
];
$LinkMap=halfmat+Transpose[halfmat];

res=N[Plus@@(Flatten[$LinkMap])/$NBATH];
Return[res];
]

SetLink[maxorder_]:=Module[{res,lkList,len,tstList,fltLKlist,fltTstList},
lkList=Table[FormLinkTable[i],{i,1,maxorder}];
tstList=Table[testLink[#,i]&/@lkList[[i]],{i,1,maxorder}];

fltLKlist=Flatten[lkList,1];
fltTstList=Flatten[tstList,1];
len=fltLKlist//Length;

Table[IsConnect[fltLKlist[[i]]]=fltTstList[[i]],{i,1,len}];
res=len;
Return[res];
]

FormLinkTable[dim_]:=Module[{res,od,digi,tri},
od=(dim-1)dim/2;
res=Table[

digi=IntegerDigits[q,2,od];
tri=Table[If[i>j,digi[[(i-2)(i-1)/2+j]],0]
,{i,1,dim},{j,1,dim}];
tri+Transpose[tri]

,{q,0,2^od-1}];

Return[res];
]

ClusterGrowth[clst_]:=Module[
{res,linklist,newlinklist,linklen,i,newelement,newClst},

linklist=Flatten[(Pick[$BathIndexList,#,1]&)/@$LinkMap[[clst]]]//Union;

newlinklist=Complement[linklist,clst];
linklen=newlinklist//Length;

res={};
For[i=1,i<=linklen,i++,
	newelement=newlinklist[[i]];
	If[isGrowthBy[clst,newelement],
		newClst=Sort[clst~Join~{newelement},Greater];
		res=Append[res,newClst];
	];
];

Return[res];
]

GenerateFullClusters[maxOrder_]:=Module[
{res,order1,i,seedclst,len,newclst,prog},

order1=List/@$BathIndexList;
FullClusters={order1};

For[i=2,i<=maxOrder,i++,
	seedclst=FullClusters[[i-1]];
	len=seedclst//Length;
	
	PrintTemporary["Generating order "<>ToString[i]<>" clusters..."];
	Monitor[
		newclst=Join@@Table[prog=j;ClusterGrowth[seedclst[[j]]],{j,1,len}];,
		ProgressIndicator[Dynamic[prog],{1,len}]
	];
	FullClusters=FullClusters~Join~{newclst};
];
FullClusters=Sort/@FullClusters;
res=Length/@FullClusters;
Return[res];
]

GenSubIdx[]:=Module[{res,fulllist,fulllength,i},
	fulllist=Join@@FullClusters;
	fulllength=Length[fulllist];
	
	For[i=1,i<=fulllength,i++,
		Pos[fulllist[[i]]]=i;
	];
	SubIndexList=subIdx/@fulllist;
	res=Length[SubIndexList//Flatten];
	Return[res];
]

isGrowthBy[clst_,newelement_]:=Module[
{order,listSmall,smallLen,i,fullList,subclst,sublinkmap,isC},
order=clst//Length;
listSmall=Select[clst,#<newelement&];
smallLen=listSmall//Length;
If[smallLen==0,
	Return[True]
];
If[smallLen==order,
	Return[False]
];

(*if 0<smallLen<order*)

fullList=Sort[clst~Join~{newelement},Greater];
For[i=1,i<=smallLen,i++,
	subclst=Delete[fullList,-i];
	(*Print["full=", fullList,"; sub=",subclst];*)
	sublinkmap=$LinkMap[[subclst,subclst]];
	isC=IsConnect[sublinkmap];
	If[isC,Return[False]];
];

Return[True];
]

testLink[lk_,nspin_]:=Module[{res,spi,seed,len,newElement,lenNew},
spi=Range[nspin];
seed=1;
newElement=Pick[spi,lk[[seed]],1];
seed={seed}~Join~newElement;


While[Length[newElement]>0,
len=Length[seed];
newElement=Complement[Union[Join@@Table[Pick[spi,lk[[seed[[i]]]],1],{i,1,len}]],seed];
seed=seed~Join~newElement;

lenNew=Length[newElement];
];

res=If[Length[seed]==nspin,True,False];
Return[res];
]

subIdx[clst_]:=Module[{res,len,subs,idxList,res1},
len=clst//Length;
If[len==1,Return[{0,{}}];];

(*len > 1*)
subs=Subsets[clst,{1,len-1}];
idxList=Pos/@subs;
res1=DeleteCases[idxList,_Pos];
res={Length[res1],res1};

Return[res];
]

End[] (* End Private Context *)

EndPackage[]